On the Josephson effect in a Bose-Einstein condensate 
subject to a density dependent gauge potential 
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We investigate the coherent dynamics of a Bose-Einstein condensate in a double well, subject to 
a density dependent gauge potential. Further, we derive the nonlinear Josephson equations that 
allow us to understand the many-body system in terms of a classical Hamiltonian that describes 
the motion of a nonrigid pendulum with an initial angular offset. Finally we analyze the phase- 
space trajectories of the system, and describe how the self-trapping is affected by the presence of 
an interacting gauge potential. 
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Introduction. The ability to manipulate individual 
quantum particles has caused a landslide of fascinating 
discoveries that have substantially increased our under- 
standing of both the microscopic and macroscopic world. 
Of principle interest are condensed matter systems where 
it is now possible to realise Feynman's vision of quantum 
simulation: the emulation of one system of interest with 
anotherp]]. Atomic condensates formed of bosons and 
fermions have been utilised to study a plethora of ideas 
due to the experimental and theoretical control they af- 
ford. One prominent example that has been the sub- 
ject of intense study over the last few years is the abil- 
ity to create artificial gauge potentials in ensembles of 
ultracold matter The ability to study gauge theo- 
ries in this context has opened up new directions in the 
ultracold-atoms landscape such as spin-orbit coupling^, 
Hall physics Q and even relativistic effects 0, 0. The ex- 
perimental control achievable over these systems means 
that one has the ability to study single and many-body 
physics with a sophisticated level of control. 

To create these "artificial" vector potentials, one can 
for instance stir the condensate with a laser, a technique 
which has been used to realise the vortex lattice of a 
condensate Q. This has its limitation in that we can only 
create a constant magnetic field in the rotating frame. 
Alternately; one can use optical couplings that lead to 
dark state dynamics [lOj or Raman transitions fll|. Fur- 
ther, one also has the option to study many-body effects 
on a lattice. To induce gauge potentials in the discrete 
setting, laser assisted tunnelling can be used in order 
to prepare the required phases for the tunnelling am- 
plitudes between the individual sites that constitute the 
lattice[H, El. 

A common feature of both the continuum and the lat- 
tice gauge theory is that they are static; in the sense 
that they are determined by the external laser couplings, 
and are not affected by the motion of the atoms. It was 
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shown in fl4| how one can induce an effective back ac- 
tion between the gauge field and the dressed states of 
the light-matter interaction, resulting in a gauge poten- 
tial that depends on the density of the quantum gas. In 
this paper we demonstrate how this continuum interact- 
ing gauge theory can be applied to a two-site lattice, from 
which we study the coherent transport between the two 
sites via the nonlinear Josephson relations which we also 
derive and numerically solve. 

A continuum interacting gauge theory. The question 
of how to study and understand the many-body problem 
lies at the heart of any realistic attempt to construct a 
theory of interacting particles in many sub-disciplines of 
physics. It is simply the case that by studying systems 
of particles with many different, interacting degrees of 
freedom one is left in a situation that is analytically and 
numerically intractable. One methodology to tackle this 
is to partition the particles according the types of motion 
that occur within the system. The most prominent ex- 
ample is given in atomic systems where the motion of the 
nucleus and the electrons are divided into slow and fast 
degrees of freedom respectively, the Born-Oppenheimer 
approximation. 

This dichotomy can be extended to other systems 
whose degrees of freedom can be separated in a similar 
fashion. Cold atom systems represent a highly control- 
lable scenario to study many-body physics at ultracold 
temperatures. However, as these atomic gases are charge 
neutral, one does not instantly have access to the many 
paradigmatic effects that pepper the condensed matter 
physics of charged particles. To redress this, one must 
find a way to simulate the mathematical structure of a 
gauge theory with these charge neutral systems. 

For bosonic systems, a simple model one can utilize is 
a system of N particles whose state space is spanned by 
the atomic states {|1), |2)}. These states along with the 
adiabatic theorem give us the ingredients to construct 
our Born-Oppenheimer like approximation. 

Our goal then is to derive an equation of motion for 
the N two-level systems at the mean-field level. As a first 
step we define the Hartree wave function for our system 
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as |\&) = <E)^Li\^i), where is the single particle wave 
function. Consequentially, one can write down an energy 
functional that will be the workhorse of the calculation, 
which is given by E = with 



2m 



1 + H lm + V + V, 



and 



fr - m 

Aim — "Tr" 







(1) 



(2) 



is the Hamiltonian for the light-matter interaction. 
Equation ([1]) above also contains the important single- 
particle external potential V that we will define later. 
Bose-Einstein condensates have particle densities that 
are typically of the order 10 13 -10 15 cm -3 , and as such it is 
appropriate to model the scattering to leading order by 
two-body zero range interactions. Hence the interaction 
matrix V in Eq. (fTJ) is given by V = (l/2)diag[<7npi + 
5i2P2,ff22/02 + <7i2Pi] and we define the population density 
in the state i as pi = \^ i\ 2 . To construct an interacting 
gauge theory, we make use of the dilute property of the 
gas to construct a perturbation theory using the atomic 
dressed states of the light-matter coupling Hamiltonian, 
which we denote |x^ 0) ) = (|1) ± e i(t> \2))/V2. We wish to 

diagonalize U + V then by treating the particle interac- 
tions as weak compared to the light-matter coupling, so 
that the chemical potential satisfies /x(r) <C Ml. Con- 
sequently, the perturbed dressed states can be written 

as 



l* ± > = lxi°>> ± ^f^ ± lx<°>> 



(3) 



and the perturbed spatially varying eigenvalues gp± ± 
Ml/ 2 now contain a contribution from the local chemi- 
cal potential of the gas, the effective scattering param- 
eter being g = (gn + 922 + 2px2)/4- We transform 
the interaction term V in Eq. |T]) into the ± basis by 
the unitary transformation U'VU, where the transfor- 
mation between the atomic and dressed basis is given by 

*le{i,2} = Ei={+,-}(^IXi 0) )*i- The two-body interac- 
tion matrix V± then reads 



9 |(5ii-522) 
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P±- 



(4) 



To build an interacting gauge theory, we define a state 
vector comprised of the two basis functions |£) = 
J2i=+ - ^i\Xi)- By projecting onto one of these states 
we assume the adiabatic theorem is valid, which requires 
that the un-projected state have negligible population. 
Thus, the effective Hamiltonian Eq. ([T]) becomes 



H+ = 



(p-A ± [r;p ± (r,t)]y 
2m 



V±(r)+E 



§P±, (5) 



with Eq = W ± Ml/2. The density-dependent geometric 
phase A-j- = ifi(x±\^X±) arises from the spatial depen- 
dence of the perturbed dressed states, and is also known 



as the mead-berry connection[15|. The scalar geometric 

phase is defined as W = ^|(x_|Vx+)| 2 . Using the def- 
inition of \x±) the vector geometric phase is then given 
to leading order by A± = A' ' ± aip±(r). There is a 
single as well as a many-body contribution to A± , where 
the single particle vector potential is A^ ) = — -|V</> and 
ai = V0(<7n — g22)/8fl determines the strength of the 
density-dependent vector potential. 

To study the dynamics of the condensate, we can derive 
a Gross-Pitaevskii like equation of motion by minimising 
the energy functional £ = ($\(iHd t - H±)\^). Without 
loss of generality we minimize with respect to "f^, and 
drop the ± subscripts on p±, ty± and A±, thus the mean- 
field equation of motion reads 



(P-A) 
2m 



ai -3 + V{r) + W + gp 



*• (6) 



Interestingly, we now have two distinct types of nonlin- 
earity appearing in Eq.©, the standard 1^1 2 term as well 
as a current j that appears at the mean-field level, given 
by 



n 

2mi 



* V 



-A )** - 

h 



V - -A I* 

h 



(7) 



As with the single particle case, the continuity equation 
that connects the probability density to the probability 
current is given by dtp + V • j = 0, although it is stressed 
that the current in Eq. ((BJ is a collective effect. The 
mean-field scalar potential W is given to leading order 
by W = |A(°) | 2 /2m. This model has recently been stud- 
ied in a one-dimensional context, where novel effects like 
chirality[3, EH are present due to the current j and the 
density-dependent gauge field in Eq. (j6]). Experimental 
realization would require atoms with long lived excited 
states; for example the transition 1 Sq Pi in Sr might 
be suitablefnj. 

One- dimensional model. Some of the most striking 
effects of coherent matter in the ultracold temperature 
regime have been elucidated with bosonic atomic conden- 
sates. Early experiments and theoretical work focussed 
on understanding interference with matter waves [l8i[l9| . 
Bragg scattering [23| and applications to matter wave 
optics such as the realization of nonlinear effects like 
solitons[2ll|. One of the most surprising yet paradigmatic 
effects in quantum mechanics is the quantum tunnelling 
of particles. For macroscopic systems this is encompassed 
by the Josephson effect: the tunnelling current that is 
produced by placing an insulating barrier between two 
particle reservoirs. The inherent nonlinearities that are 
present in the study of interacting ultracold bosonic gases 
make these systems particularly interesting, and has lead 
to the prediction and realization of effects such as self- 
trapping [23423 ■ A detailed overview of these effects is 
given by Gati et. al. [26[. Extensions of the Josephson 
effect to systems incorporating non-abelian gauge fields 
have recently been described [23, HH]. Here, we show how 
the interacting gauge theory described in the previous 



3 



section can be placed onto a two-site lattice, and study 
the population dynamics of the resulting lattice gauge 
theory (LGT). 

We envisage the situation where the cloud of atoms 
is confined so that any transversal dynamics are effec- 
tively frozen out, meaning a one-dimensional mean-field 
description of the BEC is justified. We also define <j> = kx 
as the phase of the incident laser, and also re-define the 
state <3>(x) = e~' lkx ^ 2 ip(x), which yields: 



ihdtip 



1 

2m 



(p-a 1 p) 2 + a 1 j(x) + V(x) + W + gp 



1>, (8) 



with W = H 2 k 2 /8m, and a\ = k(gu — g22)/8flS t gives 
the strength of the current nonlinearity. The parameter 
St defines the transversal area of the ID cloud. Next, we 
perform a gauge transformation on Eq. ([SJ , 



ip(x, t) = exp 



ICL\ 



dx'p(x', t) - iWt/h^(x, t), (9) 



the resulting equation of motion for &(x,t) is given by 

h 2 



ihd f <& 



-^d 2 + V(x)-2 aiJ (x)+g\$\< 



$, (10) 



where the gauge transformed current appearing now in 
Eq. (JTUJl is given by j(x) = (h/m)Im($*(x)d x $(x)). We 
can then write the one-dimensional mean-field Hamilto- 
nian in Eq. (|10jl in second-quantized form 



H = J dx &(x) ( - ^d 2 x + Vix)^ 

-2ai / dx &(x)J(x)$(x), (11) 



where the normal ordered current operator J(x) that ap- 
pears in Eq. (|11[) is given in second quantized form by 



^ ^ 2mi 



&(x)d x $(x)-d x &(x)$(x) 



(12) 



To proceed, we require a model potential V(x). Typically 
one is interested in situations where the tight binding ap- 
proximation can be made, such that the height of the lat- 
tice is greater than the chemical potential at each individ- 
ual well. Experimentally, an extended one-dimensional 
lattice can be realized by counter-propagating laser 
beams and the fact that the energy of an individual 
atom is shifted by an amount AE = —^a'(uj)(£(r,t) 2 ) t , 
where a' (to) is the real part of the dynamical polariz- 
ability of the atom, and (£(r,t) 2 ) t is the time averaged 
square of the electric field. However, as we are inter- 
ested in a two-site system, schemes involving electrostatic 



interactions [29] or atom chips [30j are more appropriate. 
Here, we consider the model potential (3lj| 



V(x) 



2b 



(13) 



The double-well potential defined by Eq. (fTB"]) above has 
its minima situated at x m i n = ±\Jmuj 2 /Ab, and close 
to these points V(x) is harmonic, an approximation we 
will use to perform the tight binding calculation. Hence 
the normalized local ground states of the left and the 
right well are given by T]i/ r (x) = (2/tto- 2 ) 1 / 4 exp(— [(x ± 

^min)/c] 2 ) respectively, and we define a = yj2h/muj as 
the width of the ground states. 

Mean-field tight binding calculation. We are now in 
a position to derive a mean-field tight binding Hamil- 
tonian describing the population dynamics between the 
two wells. To do this we use the two- mode approxima- 
tion, which entails expanding the field operator $ as 



<l = r)i(x)ci + r] r (x)c r 



(14) 



where the operator q and c r destroy particles in the left 
and right wells respectively. Note that we are assuming 
that there is a large separation between the ground state 
and the first excited state of each individual well so that 
the dynamics are well described by assuming the particles 
are located in one of the two ground states. If we insert 
Eq. (|14j) into the Hamiltonian given previously by (|11| . 
we arrive at: 



ijkl 



(15) 



ijkl 



the sums being taken over both the left and right wells. 
The parameters Jy, Uijki and Xijki are overlap integrals, 
and give us a way to identify the most important terms in 
Eq. (TT5|) . The overlap integrals can be readily calculated 
using our definition of the ground states r)i/ r (x) given 
previously, thus one finds that the lattice Hamiltonian 
Eq. [T51 simplifies to 

H = — J{c\c r + c\ci) + U{hi(hi - 1) + h r (h r - 1)) 
+ Ti(c\jci + cljcr) + T 2 (cjjc r + cljci), (16) 

where fii = cjc» is the number operator for site i and 
the constants that appear in Eq. (|16[) are given by 

U = g/V^v, 2 ri = 8^ix min e- 3a; «"n/ ff2 /^ 3 ; T 2 - 
Sha\x m x n e~ 4x ^' a I y/na 3 and the discrete current op- 
erator for the two site lattice we define by j = — i(c}c/ — 
cjc r ). Equation (|16p gives us a way to understand the 
effect of population dynamics between the two sites. It 
comprises the usual on-site interactions that appear in 
the Bose-Hubbard model given by the terms proportional 
to U , as well as the unconventional terms proportional 
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(b) 



2 2.5 3 



(c) 



(d) 



FIG. 1: (Color online) Numerical solutions to Eq. (|17|l and 
(|18[) . The population |q| 2 , ]c r ] 2 and |ci| 2 are given by 
the blue dashed line, solid red line and dashed black lines 
respectively. Figure (a) and (c) show the population oscilla- 
tions for 2U/J = 10 and Fi/J — 1, the insets in each figure 
show the current J(t) as a function of time. Both (a) and 
(b) have the initial conditions C;/ r (t = 0) = 1/V2. Figure (b) 
and (d) are plotted for 2(7/ J = 1 and Ti/J = 5. Figures (c) 
and (d) have the initial conditions ci(t = 0) = e iv ^ 2 /y/2 and 
c r (t = 0) = \jsf2. The units of time are h/J. 



to Ti which originate from the current operator, Eq. 
(fT2")l . To study the population oscillations between the 
two wells we work with the operators ti in the Heisen- 
berg picture. The equations of motion are then given by 
ih£i = [ti, H], which gives 



., dci 
in—— 
dt 



ih- 



dc r 
~dt 



— Jc r + lUhiii + Tx{jci + i{hi + n r )c r ) 

+ r 2 (jc r + icjc r c r + ih r ci), (17) 

— Jci + 2Un r c r + Ti(jc r — i(ni + n r )ci) 

+ r 2 (jc; - i&lcidi - ihidr). (18) 



To gain an understanding of Eq. (fTT)) and (fT5)) we assume 
that the number of particles in both wells is so large that 
the operators q and c r may be treated as classical quanti- 
ties, and as such can be replaced by their expectation val- 
ues Q and c r . Figure [T] shows the numerical solutions to 
these equations in different parameter regimes. In Figure 
1 (a) and (c) the parameters 2U/J = 10 and T±/J = 1 
were used, whilst for (b) and (d) 2U / J = 1 and T\/J = 5. 
For figures (a) and (b) the initial phase different was 
di - 9 r = A6 = 0, but for figure (c) and (d) A6 = tt/2 
was used. We further assumed that x m i n /a 3> 1, so that 
r 2 <C Ti. Figure Q] shows the Rabi like oscillations that 
are synonymous with two state systems. We see in (b) 
the increased current strength has caused the speed of the 
population oscillations to increase. On the other hand, 
Fig. [1] (c) shows how strong on-site interactions changes 



the dynamics, the current showing an unusual 'dip' (see 
inset). Finally Fig. Q] (d) shows how the dynamics are 
reduced when there is an initial phase difference of 7t/2 
between the sites and the current strength is stronger 
than the Hubbard term. 

Phase-space analysis. Figure [T] shows how the current 
nonlinearity in Eq. (JT5J) affects the population oscilla- 
tions between the two wells. To investigate the properties 
of this system further, we can re-cast the variables of the 
problem in terms of the population difference between the 
two sites and the phase difference. This methodology has 
previously been utilised to show how a charge neutral in- 
teracting BEC in a symmetric double well potential can 
be understood in terms of a nonrigid pendulum^, HH, 
and gives an intuitive way to study the phase-space prop- 
erties of the many-body system[25|. We begin by assum- 
ing again that the number of particles is so large that 
the operators that appear in Eq. (TTrl)) can be substi- 
tuted for their respective eigenvalues, which in turn can 
be replaced by the polar variables Cj = \fN~ie %6i , where 
i = l.r. If we define the population and phase differences 
by z(t) = (Ni - N r )/N t and ip(t) = 6 r - 6 t respectivly, 
then we obtain a classical Hamiltonian from Eq. (|16l) 

Az 



H =— y/l — z 2 cos(tp) — 7i\A — z 2 sm((p) 



(19) 

where the dimensionless parameters read A = UN 2 / J, 
ji = T t N 2 /J and AE = UN 2 /2J. The tunnelling cur- 
rent is j = —Nt\/l — z 2 sin(<p). The variables z and ip 
are canonically conjugate, with z = and tp = ^j, 

which gives the coupled equations of motion 



z = — \fl — z 2 sm(ip) + 7i yl — z 2 cos(ip) 

+ 2 72 (l-z 2 )cos(2^), (20) 



: cos(ip) + Az + 7i 



sin(<^) 



■ 27 2 zsin(2^?). 



(21) 



The Hamiltonian Eq. (|19|) and the nonlinear Joseph- 
son equations Eq. (|2U|) and (f2"l"T) allow us to understand 
the dynamics and phase-space properties of the two site 
model. In particular, we note that Eq. can be un- 
derstood in terms of a classical nonrigid pendulum. Our 
model differs from that presented in [22j by the addi- 
tional term proportional to 72. Hence, if we drop the 72 
term from Eq. ([191) we can map the classical Hamilto- 
nian onto the nonrigid pendulum model as presented in 
[22l |. The result is that the phase angle (p has an extra 
initial offset term due to the current nonlinearity of the 
underlying continuum model. In this limit Eq. (|19[) can 
be simplified to 



H 



Az 2 n r- 



z 2 cos(<^ -tpo), 



(22) 



where R = yl + Ti an< ^ ^he angle ip = arctan^), and 
we set AE = without loss of generality. Figure [2] shows 
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FIG. 2: (Color online) Lines of constant energy, with z vs. 
if. The initial conditions were z(Q) = 0.5 and (p(0) = it. The 
figure and its inset demonstrate two regimes: the main figure 
has A = 2 and 71 = {0,0.5,2}. The inset has A = and 
71 = {0, 1, 2}, where the different values of 71 correspond to 
the smallest through to the largest curves in the figures above. 

the phase-space trajectories of variables z, cp. In Fig. [2] 
we are interested in two limits, the first being the phase- 
space with on-site interactions, A = 2, whilst the inset 
is plotted for no on-site interactions, A = 0. This is 
justifiable as we could for example use Feshbach reso- 
nances in order to achieve g = 0. The plots for A = 2 
show how increasing the strength of the current causes 
the curves obtained from the numerical solutions to equa- 
tions ([20]) and (|21|l to increase in size. The inset shows 



how increasing the strength of the current without on- 
site interactions gives a displacement of the curves by an 
amount ipg. To further quantify this, we calculate the 
critical value A c that determines the point at which the 
nonrigid pendulum is given an initial kick that pushes 
it over the vertical <p = n point. This is given by the 
condition H[z(0),<p(0)] > 1 so that 

_ 2(1 + y/l + z(0)*Rcos[tp(0)-tp ]) 

If we choose A = A c we are in the situation whereby the 
pendulum is in an upright position, so that the period of 
oscillation diverges. 

Conclusions. It was shown how a one-dimensional 
mean-field theory describing density dependent gauge 
fields can be cast into a two-site tight binding model, 
using a symmetric double-well potential. It was seen 
that this discrete formulation also features a current op- 
erator as well as the usual on-site interactions that ap- 
pear in the standard Bose-Hubbard model. Finally, a 
phase-space analysis was presented, by way of the non- 
linear Josephson equations. It was found that a classi- 
cal Hamiltonian can be written down that describes the 
motion of a nonrigid pendulum, with an initial angular 
offset that depends on the strength of the density depen- 
dent gauge potential. Further analysis of this model is 
warranted. For instance, the unconventional transport 
mechanisms could have a significant impact on the su- 
perfluid behaviour, and is likely to affect non-trivially 
the elementary excitations of the system. 
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